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I. INTRODUCTION 

The variational and diffusion quantum Monte Carlo 
methods (VMC and DMC) are zero-temperature stochas- 
tic techniques for evaluating the expectation values of 
time- independent operators [ll-Q- These methods are 
particularly well suited for calculating the ground-state 
energies of large assemblies of interacting quantum par- 
ticles. The central object is an approximate trial wave 
function whose accuracy controls the computed energy 
and the intrinsic statistical fluctuations in the calcula- 
tions. It is therefore of particular importance to develop 
accurate trial wave functions. 

Expectation values in VMC are evaluated using 
importance-sampled Monte Carlo integration. In DMC 
the ground state is projected out by evolving the 
Schrodinger equation in imaginary time. Such projector 
methods suffer from a fermion sign problem, in which 
errors in the propagation increase exponentially in imag- 
inary time as the algorithm amplifies any spurious com- 
ponent of the lower-energy bosonic state. This problem 
may be evaded in DMC by employing the fixed-node ap- 
proximation in which the nodal surface is fixed to 
be that of a suitable trial wave function. The resulting 
DMC energy is greater than or equal to the exact energy 
and less than or equal to the VMC energy computed with 
the same trial wave function. The VMC energy depends 
on the entire trial wave function, but the DMC energy 
depends only on the nodal surface of the trial wave func- 
tion. 

One of the appealing features of VMC and DMC is 
that virtually any form of trial wave function can be 
used. The main criteria are that the wave function must 
obey the correct symmetry under particle exchange, it 
should be flexible enough to describe the system of inter- 
est, and that it should be possible to evaluate it rapidly. 
The analytic properties and normalizability of the trial 
wave function must be such that the energy expectation 
value is well-defined. The simplest fermionic wave func- 



tion is a Slater determinant, which describes exchange 
but not correlation. Multidetcrminant wave functions, 
pairing wave functions such as geminals [H , and backfiow 
transformations can also be used. The most fruitful 
method of going beyond the Slater determinant is, how- 
ever, to multiply it by a Jastrow factor which leads 
to the Slater- Jastrow wave function. The Jastrow factor 
is normally chosen to depend on the interparticle sepa- 
rations, which introduces correlation into the wave func- 
tion. The introduction of a Jastrow factor often leads to 
the recovery of 80% or more of the correlation energy of 
electronic systems 

The Jastrow factor is chosen to be everywhere positive 
and symmetric with respect to the exchange of identical 
particles in order to maintain the nodal surface defined 
by the rest of the wave function. One of the features 
of the Jastrow factor is that it can conveniently be used 
to enforce the Kato cusp conditions Q , which determine 
the behavior of the wave function when two charged par- 
ticles approach one another. Enforcing the Kato cusp 
conditions does not necessarily improve the variational 
energy, but the reduction in the statistical fluctuations 
in the energy is often very important. 

DMC can be viewed as VMC with a perfect Jastrow 
factor, but improving the Jastrow factor can improve 
DMC calculations in several ways. The DMC algorithm 
is subject to time-step errors and to (normally very small) 
population-control errors @ that are reduced by improv- 
ing the trial wave function. Evaluating expectation val- 
ues of operators that do not commute with the Hamil- 
tonian is not straightforward in DMC, but using highly 
accurate trial wave functions helps in achieving more ac- 
curate results. Similar considerations apply when us- 
ing nonlocal pseudopotentials, which involves making ap- 
proximations that are ameliorated by improving the trial 
wave function [lol . [ll| . As the fundamental limitation on 
the accuracy of DMC is the quality of the nodal surface, 
it is desirable to use trial wave functions with optimizable 
nodal surfaces as afforded by, for example, multidetcrmi- 
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nant wave functions and backflow transformations. A 
good Jastrow factor can account for the bulk of the dy- 
namical correlation energy, which allows the optimization 
of parameters that affect the nodal surface to achieve a 
better nodal surface. 

Here we introduce a highly flexible form of Jastrow 
factor which allows for the introduction of a variety of 
terms involving arbitrary numbers of particles. Our main 
motivation is to be able to implement quickly different 
functional forms and explore the importance of different 
correlations in any physical system we study. 

Jastrow factors correlating several electrons have been 
used in earlier calculations, such as those of Refs. 
[T^ . We study the effects of various three-body Jas- 
trow terms and introduce a four-body van der Waals- 
likc term. We also construct anisotropic Jastrow factors 
that can capture the natural symmetries of a system. We 
have successfully applied these Jastrow factors to a va- 
riety of systems, and we report results for the one- and 
two-dimensional homogeneous electron gases, the Be, B, 
and O atoms, and the BeH, N2, H2O, and H2 molecules. 
Our VMC and DMC calculations were performed using 
the CASINO package [l^. Hartree atomic units are used 
throughout (h = |e| = me = 47reo = 1). The structure of 
the paper is as follows. We describe the form and prop- 
erties of the general Jastrow factor in Sec. HT] Specific 
examples of the construction of Jastrow terms are given 
in Sec, mil and results obtained using them are presented 
in Sec. IIVI Finally, we draw our conclusions in Sec. IVl 
Implementation details are given in Appendix \^ and an 



Slater determinants, with different advantages and disad- 
vantages depending on the system. 

Typically J(R) is constructed as a sum of terms, e.g.. 



example can be found in Appendix 
summaries of our data in this paper 



We report only 



II. CONSTRUCTION OF A GENERIC 
JASTROW FACTOR 

Quantum Monte Carlo (QMC) methods can be ap- 
plied to systems which can be generically described as 
an ensemble of N quantum particles and M sources of 
external potential. The most common type of QMC sim- 
ulations are electronic calculations, where the quantum 
particles are electrons and the sources of external poten- 
tial are fixed nuclei (or pscudopotentials) . For simplicity, 
we refer to quantum particles as electrons and to external 
potentials as nuclei in the rest of this paper. Our Jastrow 
factor is applicable to other types of quantum particles 
and external potentials. 

Any trial wave function can be written in the form 



4't(R) = e"'(^)*s(R) , 



(1) 



where ^'s(R) is the part of the wave function that im- 
poses the symmetry and boundary conditions, and e''^^^ 
is a Jastrow correlation factor which is constrained so 
that the symmetry and boundary properties of \I's(R) 
are transferred unmodified to 4't(R)- In the single- 
determinant (SD) Slater- Jastrow wave function, >I's(R) 
is a Slater determinant. There are various alternatives to 



J(R) = Jc-c(R) + Jo-n(R) + Jc-o-n(R) 



(2) 



where "e-e" stands for "electron-electron," "e-n" for 
"electron-nucleus," etc. Each of these terms involves dif- 
ferent numbers of electrons n and nuclei m. We shall refer 
to n and m as the electronic and nuclear ranks of a term, 
respectively, which are constrained to satisfy n + m > 2^ 
n > 1, and to > 0. We have designed a generic Jas- 
trow term of selectable ranks, J„^m(R), such that the 
total Jastrow factor is constructed as the exponential of 
a sum of one or more terms of the desired ranks. In this 
notation, Jc-c = Jc-n = Ji.i, etc. 

The function Jn,„i(R) is a sum over all sets of n elec- 
trons and m nuclei in the system of a parameterized func- 
tion of the e-e and e-n relative position vectors within 
each such set. While alternatives exist, a natural way 
of parameterizing this function for arbitrary values of n 
and m (implying an arbitrary number of variables in the 
function) is to expand it in products of functions of the 
individual e-e and e-n vectors. Thus, we construct our 
Jastrow factor using pairwise objects as building blocks, 
and in what follows we describe these objects and derive 
the properties of J„^m(R) that follow from those of the 
pairwise objects. 

We name the e-e functions used in the expansion 
$^(r), where r is the relevant e-e relative position vec- 
tor, v is the index of the function within a chosen ba- 
sis of functions, and P is the e-e dependency index, 
which allows the use of different optimizable parame- 
ters, if present, for parallel- and antiparallcl-spin electron 
pairs, for example. Similarly, the e-n basis functions are 
0^(r), where r is the relevant e-n relative position vec- 
tor, fi is the index of the function within the chosen basis 
set, and S is the e-n dependency index of the basis set, 
which allows the use of different parameters for up- and 
down-spin electrons around a given nucleus, or for dif- 
ferent atoms, for example. In the case of nonelectronic 
systems, e-e and e-n dependency indices arc used to dis- 
tinguish between particle types and spins. 

We introduce a compact notation for defining J„.m(R). 
We represent the n electronic indices by the integer vec- 
tor i = {ii, 12, ■ ■ . , in}, each of whose components takes 
a distinct value between 1 and N , and the to nuclear in- 
dices by the integer vector I = {Ji, /2, . . . , /,„}, each of 
whose components takes a distinct value between 1 and 
M . For each term in the Jastrow factor wc define the e-e 
and e-n dependency matrices P and S of respective sizes 
N X N and x M containing the dependency indices 
Pij and Sii for each e-e and e-n pair. The components 
of P and S can be made equal depending on the symme- 
tries of the system, including particle distinguishability 
and geometrical symmetries which make different nuclei 
equivalent. 

Likewise, it is convenient to use matrices to represent 
the basis functions involved in the Jastrow factor term. 
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For c-e basis functions, each row and column of the nxn matrix $ corresponds to an electron, 
I 



$^(i) 







(3) 



We refer to the nxn matrix formed by the e-c depen- 
dency indices {^i„i^}a,^i=i,...,n asP(i). BothP(i) and the 
nxn matrix of e-e expansion indices v are defined to be 
symmetric, and this fact has been used in Eq. ([3]). Noting 
that Tji = —Vij, and restricting the e-e functions to be 
either symmetric or antisymmetric about the origin, one 
finds in Eq. that matrix $ is symmetric, antisymmet- 
ric, or asymmetric depending on whether the functions 
in the basis set are all symmetric, all antisymmetric, or 
both types are present, respectively. 

For e-n basis functions each row of the nxm matrix 
corresponds to an electron and each column to a nucleus. 



0|(i,I) 



(4) 



We refer to the nxm matrix formed by the e-n depen- 
dency indices {^i^/^ /3=i,....m as S(i, I), and the 
nxm matrix of e-n expansion indices is fj,. 

We write Jn,m as a sum of contributions from each 
group of n electrons and m nuclei in the system. 



s.v. s.v. 



Jn.m — j ^ ] y ] t/ri,m(i, I) — ^ ] ^ ] Jn,rn{i, I) 7 (5) 



where summations with vector indices represent sums in 
which every component of the vector is a summation in- 
dex, and "s.v." (for "sorted vector") indicates that the 
sum is restricted to vectors whose components are sorted, 
e.g., ii < io < . . . < in, which avoids redundant contri- 
butions |17| . The contribution from the n-electron and 
m-nucleus group {i, 1} is 

I) = E E >^'^''' n n I) ' (6) 

where A are the linear parameters, summations with ma- 
trix indices represent sums in which every component of 
the matrix is a summation index, Y[ acting on matri- 
ces implies the product of all of their components, and 
"u.t." means that the relevant operation is restricted to 
the upper-triangular portion of the e-e matrices involved, 
excluding the diagonal. 
A. Symmetry properties of the linear parameters 

Equation ^ imposes the condition that J„ „,(i,I) 
must not depend on the specific ordering of the electrons 
and nuclei listed in i and I. Let A and B be permuta- 
tion matrices of respective sizes nxn and m x m such 
that Ai and BI arc integer vectors containing reordered 
electronic and nuclear indices. The value of J„.„i(Ai,BI) 
should therefore equal that of Jn,m{h^), 



j„..(Ai,Bi). E E Agijii! n^(Ai)n0|(Ai,Bi) 

K(Ai) m(BI) 

= E E At^!;^A^i.^'''^^nAti(i)A-nAe|(i,i)B- 



\AP(i) A^,AS(i,I)B^ 



$i(i),Ajn^i(i)n4(i'i)' 



(7) 



where 



n""At£(i)A^ 

n"'*F(i) 



(8) 



which is for basis sets consisting only of symmetric 
functions, while in the presence of antisymmetric basis 
functions it may be +1 or —1 depending on the precise 
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permutation performed by A. Equating the right-hand 
sides of Eqs. ^ and (O one finds that 



\P(i).S(iJ) ^ 
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\AP(i)AT,AS(i,I)BT 



This equation represents the basic symmetry property of 
the hnear parameters of the Jastrow factor, which im- 
phes that a parameter with a given set of superindices 
{P(i), S(i, I)} is determined by another parameter with 
a permuted set of superindices {A P(i) A""", A S(i, I) B"""}. 
This redundancy is removed by considering only one of 
the possible permutations of {P(i), S(i, I)}. We call this 
particular permutation of {P(i), S(i, I)} the signature of 
the group of particles {i,I}, 

{P(i),S(i,I)} = {UP(i)U'^,US(i,I)RT} , (10) 

where the permutation matrices {U, R| are computed by 
applying a matrix-sorting algorithm [l8| to {P(i), S(i, I)}. 
In our terminology, the set of linear parameters whose 
superindices reduce to the same signature constitute a 
parameter channel. Only those parameters whose su- 
perindices equal the signature of a channel need be 
stored, and any other linear parameters in the channel 
can be computed from them via Eq. ©. 

The signature {P(i), S(i, I)} may contain repeated en- 
tries such that there exist permutation matrices {A, B} 
that leave the signature unchanged, 

{P(i), S(i, I)} - {A P(i) AT, A S(i, I) B^} , (11) 



in which case Eq. ^ becomes 



;^P(i).s(iJ) 
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$^(i),A 



xP(i),S(i,I) 



(12) 



Equation is the symmetry constraint that relates 
linear parameters within a channel, which can be imposed 
as detailed in Sec. Ill CI 



B. Indexing of basis functions 

The components of and n are the e-e and e-n expan- 
sion indices. We define the expansion indices so that they 
can each take any value between 1 and the e-e expansion 
order p, and between 1 and the e-n expansion order q, re- 
spectively. We factorize an optional cutoff function into 
and so that 



<I.r(r)=r(r)0r(r), 



for v > 0, and 



e^(r) 



/(r)^^f(r) , 



(13) 



(14) 



for ^ > 0, where and are the e-e and e-n cutoff 
functions and (j)^ and 6^ are functions from a suitable 
basis set. This factorization allows an efficient imple- 
mentation of localized Jastrow factor terms. 



Additionally, we allow expansion indices to take a value 
of zero with the special meaning that $f (r) = Qq (r) = 1 
for all P, S, and r. Note that these 0-th functions do 
not contain cutoff functions. This allows us to construct 
terms with specialized functional forms, such as those 
involving dot products of vectorial quantities. 



C. Constraints 

Constraints on the parameters can be expressed in the 
form of a system of equations involving the linear param- 
eters and the basis function parameters. We restrict our 
analysis to linear constraints on the linear parameters, 
and constraints that can be imposed on the nonlinear 
parameters contained in a basis function independently 
from the linear parameters and from nonlinear parame- 
ters in other basis functions. 

Linear constraints on the linear parameters can be 
imposed using Gaussian elimination, as described in 
Ref. 01 • The matrix of coefficients may depend on the 
nonlinear parameters in the basis functions, if present, 
and the linear system is usually underdetermined, result- 
ing in a subset of the parameters being determined by the 
values of the remaining parameters, which can be opti- 
mized directly. 

When a constraint results in setting specific linear pa- 
rameters to zero, it is more convenient simply to remove 
them from the list of linear parameters. This is accom- 
plished by disallowing the indices and ^ from taking 
the values corresponding to linear parameters. We call 
this an indexing constraint. 



1. Symmetry and antisymmetry constraints 

Symmetry constraints must always be imposed, oth- 
erwise the trial wave function is unphysical and cal- 
culations give erroneous results. Symmetry constraints 
amount to equalities between pairs of parameters as per 
Eq. (IT2]) . When two of these equalities relate the same 
pair of parameters with opposite signs, e.g., Ai = A2 and 
Ai = — A2, which implies Ai = A2 = 0, both parameters 
arc eliminated using indexing constraints. 



2. Constraints at e-e and e-n coalescence points 

The Coulomb potential energy diverges when the posi- 
tions of two electrons or an electron and a nucleus co- 
incide. However, the local energy of an eigenstate of 
the Hamiltonian, including the exact ground-state wave 
function, is finite and constant throughout configuration 
space. Divergences in the local energy are therefore not 
a feature of the exact wave function and can lead to poor 
statistics in QMC calculations; hence it is important to 
avoid them. The kinetic energy must diverge to cancel 
out the potential energy and keep the local energy finite. 
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which is achieved by demandingthat the wave function 
obeys the Kate cusp conditions Q . For any two charged 
particles i and j in a two- or three-dimensional system, 
these are 



1 d'^ 



d±l 



(15) 



where ^ denotes the spherical average oi'^ , q represents 
charge, = mirrij / {rrii + rrij) is the reduced mass, m 
represents mass, d is the dimensionality, and the positive 
sign in the denominator is for indistinguishable parti- 
cles and the negative sign is for distinguishable particles. 
Fixed nuclei are regarded as having an infinite mass. Di- 
vergent interactions other than the Coulomb potential 
would give rise to different expressions. 

It is common practice to impose the e-n cusp condi- 
tions on X&s and the e-e cusp conditions on the Jastrow 
factor. This is because typical forms of explicitly de- 
pend on the e-n distances but not on the e-c distances. 
Our implementation allows the option of applying both 
types of cusp conditions to the Jastrow factor, which 
gives flexibility in the choice of ^'s and its properties. 
In particular, we impose the cusp conditions on a sin- 
gle Jastrow factor term, and constrain all other terms 
in the Jastrow factor so that their contribution to the 
local kinetic energy is finite at e-e and e-n coalescence 
points. For nondivergent interaction potentials, such as 
most pseudopotentials, we simply require that the kinetic 
energy remains finite at coalescence points. Our imple- 
mentation is also capable of not applying any constraints 
at e-c and e-n coalescence points since this is advanta- 
geous in some cases [l^ [13] ■ 

Imposing that the kinetic energy be finite at coales- 
cence points is nontrivial if the Jastrow factor contains 
anisotropic functions. Consider the exponent of a Jas- 
trow factor J near a point where two particles coalesce, 
be it two electrons or an electron and a nucleus. The 
dependence of J on coordinates other than those of the 
coalescing particles should be smooth in the vicinity of 
the coalescence point, and therefore one should be able 
to write J = J(r), where r is the difference between the 
position vector of the two particles, and all remaining 
particles are held fixed. 

The local kinetic energy is computed from two esti- 
mators, one involving VJ(r) and the other V^J(r). We 
require both quantities to remain finite as r — 5- 0. We ex- 
pand the Jastrow factor in spherical harmonics, J(r) = 
EZo EL-/ -^^'"^ with JC'™) = fi,m{r)Yi,^ie, 4>), and 
the gradient and Laplacian of jC'™) are 

vj('''"> = //,™(r)y,,™(0,</.K 



r 86 

fl,m{r) dYi^ra{e,dp) 

r sin 6 dd> 



(16) 



y2 j(;,m) 



i{i + i)hAr) 



Yi,m{e,(b). (17) 



Let us assume that fi^m{f) is finite at the origin and ex- 
pand it to second order about r = 0, fi.rn{r) ~ /i.m(O) + 
//,m(0)^ + //'rn(0)^V2. Substituting into Eqs.'dlll) and 
([T7|. and ignoring contributions of 0{r) or higher, we 
arrive at 



+ 



r 

fi.r.m 

r sin 9 



dY, 



80 
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dY, 



d(j) 



-Ug 



3-^).//.',„(0) 



//.™(0) 



^^iWo) 



y,,„(^,(^). (19) 



The coefficient of the negative powers of r in Eqs. (|T8)) 
and ([T^ must vanish for V J and J to be finite at the 
coalescence point. This gives rise to two conditions: (a) if 
I ^ then /i,,„(0) = 0; and (b) iil^l then //,„(0) = 0. 

Application of the Kato cusp or finite kinetic energy 
constraints requires the construction of a linear system 
for each linear-parameter channel in each term of the Jas- 
trow factor based on the above equations. Let Vi^m be an 

operator such that Vi^m YaLo J2L=-i fi..m{r)Yi^rn{Q, 4>) = 
fi,m{'r)- The cusp equations associated with the coales- 
cence of electrons i and j have the form 



u.t.jO.c. n.c. 



0,0 



dri. 



(20) 



nj=0 



where Tij is the right-hand side of Eq. p^ . The label 
"o.c." (for "one contribution") denotes that the sum is 
restricted to values of the v_ sets such that the elements 
of the upper triangular portion of $^(i) are all 1 except 
that corresponding to the electron pair formed by i and 
j, and the label "n.c." (for "no contribution") denotes 
that the sum is restricted to values of ^ such that all 

elements of 0^(i, I) are 1. These restrictions are trivially 
satisfied by e-e terms. 

Parameters that do not contribute to Eq. (pO)) should 
be set by the condition that the kinetic energy docs not 
diverge at e-e coalescence points, resulting in 



u.t.,e.p. c.p 



EAf:; 



,=0 



(21) 
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for / ^ 0, and 



u.t.,c.p. c.p. 



P,S 



l.ra 



, (22) 



for I ^ \. The anisotropy of ^ul] (r^ ) at r^j = deter- 



mines which of Eqs. (PT|) and Eqs. (22) need be imposed. 
The label e.p. (for "equal-product" ) denotes that the sum 
is only over indices associated with electrons i and/or j, 
and these indices take only values such that the prod- 
uct of the pair of functions associated with Vi}^ and Vj}^ 
{flu and iiji in the e-n case) is equal throughout the sum. 
Each set of distinct two-function products and each value 
of {l,m) correspond to different equations, and each set 
of possible values of the indices not summed over corre- 
sponds to a separate set of equations. 

For the coalescence of electron i and nucleus / the cusp 
conditions take the form 



u.t.,n.c. o.c. 

E E^f. 



dVLm [e^:j(r./)] 



r^i , (23) 



while the kinetic energy is kept finite if 

U.t.jC.p. 



E %l\T'i^rn\e^(ru)]l^^^, = 0, (24) 



for Z 7^ and 



u.t.,e.p. 

E 



dru 



= 



(25) 



for 1^1. The equal-product constraint on the sum is 
now such that the sum is only over e-e indices asso- 
ciated with electron i and e-n indices associated with 
nucleus /, and these indices only take values such that 
the product of the pair of functions associated with 
and fiki is equal throughout the sum. Again, each set of 
distinct two-function products corresponds to a different 
equation, and each set of possible values of the indices 
not being summed over corresponds to a separate set of 
equations. 

Note that the equal-product constraints in the sums 
of Eqs. (PO)) and ([^ reduce to the equal-product con- 
straints described for the e-e-n / term of the Drummond- 
Towler-Needs (DTN) Jastrow factor in the Appendix of 
Ref. when natural-power basis functions are chosen. 



3. Other constraints 

It is possible to construct terms containing dot prod- 
ucts by using appropriate constraints. For example, con- 
sider the basis functions 6i(r) = x, 82(1") = y, and 
03 (r) = z. In an e-n-n term we can restrict the indices 



so that ji takes only the values (1 1), (2 2), (3 3), so that 
the contribution of electron i and nuclei / and J is Yu-Yij, 
provided we also apply a linear constraint that equates 
the three nonzero linear coefficients. Section fill Bl gives 
a practical example of a term containing dot products, 
which is used in Sec. IIV C 41 

It is also possible to introduce Boys-Handy-style in- 
dexing, |2l| where the sum of all e-e and e-n indices is 



restricted to be less than or equal to some fixed integer 
I. This is accomplished by setting the e-e and e-n ex- 
pansion orders to I and then eliminating the parameters 
that violate the conditions via indexing restrictions. 



III. BASIS FUNCTIONS AND TERMS 



Basis sets and cutoff functions 



Possibly the simplest basis set is the natural powers. 



(26) 



as used in the DTN Jastrow factor for the localized m, 
X, and / terms 0]. These functions need to be cut off 
at some radius L, for which purpose the DTN Jastrow 
factor uses the polynomial cutoff function 



D{v) = (r - Lfe{L - r) , 



(27) 



where L is an optimizable parameter, C is a positive 
integer, and 0(r) is the Heaviside step function. We also 
use a slightly different version of this cutoff function. 



P(r) = {l-r/Lfe{L-r) 



(28) 



which should be numerically superior to D{r). 

A particular variant of P{r) is the anisotropic cutoff 
function 



Air)^il-r/LfeiL-r)J2c,Y[ 



(29) 



where L is an optimizable parameter, C is a positive in- 
teger, d is the dimensionality of the system, {up}i3=i^...^d 
are unit vectors along d orthogonal directions, {ci} are 
real-valued constants, and pjj'^ are integer exponents, 

which arc constrained so that p^*'' is the same for 
all values of i. This cutoff function is simply the product 
of an isotropic cutoff function and a spherical harmonic. 
For example, with ci = 3, C2 = —1, p^^^ = (2 1 1), 
and p(^^ = (031), and the vectors pointing along the 
Cartesian axes, we obtain 



A{r) = (l-r/LfeiL-r) 



{3x^ -y^)yz 



(30) 



which is proportional to a real spherical harmonic with 
/ = 4. The advantage of describing anisotropy in the cut- 
off function rather than in the basis functions is that the 
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common spherical harmonic can be factorized out of the 
sum over expansion indices, which reduces the computa- 
tional cost. We allow different orientations to be used for 
different e-e or e-n dependency indices, which is useful 
to adapt the functional form to, e.g., the geometry of a 
molecule. 

An alternative to the natural-power basis in finite sys- 
tems is a basis of powers of fractions which tend to a 
constant as r — > cxj, and therefore do not need to be cut 
off. We define the basis 



Mr) 



(31) 



where a and b are real-valued optimizablc parameters. 
Similar basis sets with 6=1 have been used in the liter- 
ature, often in conjunction with Boys-Handy-style index- 
ing [l^ [2l| - |23j , and this basis was used in Ref . [13] with 
an early implementation of the Jastrow factor presented 
here. 

In extended systems it is important to use a basis that 
is consistent with the geometry of the simulation cell and 
has the periodicity of the system, such as a cosine basis. 



C,(r)= J2 cos(G-r) 

GGi^-th star 



(32) 



where the G vectors are arranged in stars defined by the 
cell geometry. This basis is used in the DTN Jastrow 
factor for the extended p and q terms. 

A suitable basis set for building specialized terms con- 
taining dot products is 

V^{r) = ^INTf(t^-l)/d] • UMOD(^-l,d) + l ^ ^gg^ 



than one electron and one or more nuclei are involved, we 
choose not to apply e-e cutoff functions, relying instead 
on the e-n cutoffs to fulfill this role. Additional 7V„,m 
terms used in this paper that were not part of the DTN 
Jastrow factor are A^i,2, -^3,0, A^i,3, -^2,2, A^a.ij ^^'^ ^4,0- 
In Nn,m we typically use a truncation order in the cutoff 
function of C = 3. 

We use A^j'j^ to refer to the anisotropic variant of 
Nn.m- The m term consists of natural power ba- 
sis functions N^, and the anisotropic cutoff function A, 
and "s.h." is a placeholder for describing the spherical 
harmonic. For example, for the highly anisotropic N2 
molecule we use terms such as Af ^, Af A2 i, and A2 i ■ 

In finite systems we also use the F^, basis functions 
in terms without explicit cutoff functions which we call 
Pn,m, or whcu wc forcc & = 1 in the basis functions. 
In some systems it is useful to apply Boys-Handy-style 
indexing to i^^^j, and we refer to the resulting term as 

In extended systems wc make use of the cosine basis 
functions Cjy in terms denoted Cn,rm where we choose 
expansion orders so that at least as many G vectors as 
electrons in each spin channel are included in the expan- 
sion. 

To test the fiexibility of our implementation we have 
designed an e-e-n-n Jastrow term for describing the 
correlations associated with van der Waals interactions, 
which we call V2,2- This term is capable of distinguishing 
between configurations where the electron-nucleus rela- 
tive position vectors r^/ and Tjj are parallel from those 
where they are antiparallel. Introducing a dot product 
achieves this effect, and V2,2 has the following functional 
form. 



where d is the dimensionality of the system and 
{up} j3=i....^d arc unit vectors parallel to the d Carte- 
sian axes. A term constructed using these functions 
with appropriate index-restriction constraints would con- 
sist of dot products between two vectors multiplied by a 
natural-power expansion in their moduli. 



B. Terms and notation 

We employ a condensed notation to refer to Jastrow 
terms that use certain basis functions, cutoff functions 
and constraints. Each term is represented by a single 
capital letter, with n and m as subindices. Any other 
relevant information is given as a superindex. Typically 
we use expansion orders p and q of 7-9 for two-body 
terms, 4-5 for three-body terms, and 2-3 for four-body 
terms, except when indicated otherwise. 

For simple Jastrow terms we use the natural power ba- 
sis functions and the polynomial cutoff functions P 
or D. We refer to these terms as Nn.m- -^2,0: -'^^i,!: and 
iV2,i are the equivalent of the DTN u, x, and / terms, re- 
spectively. In the iV2.i term, and in any term where more 



N Ai p q 

xA^.., {nj)N^^, iru)N^,j {rjj)ru ■ v,., . (34) 

We require basis functions to be scalars in our Jastrow 
factor, so the dot product is separated into its compo- 
nents. Hence, wc construct the T^2,2 term using Vu for 
the e-n basis with P as the e-n cutoff functions, and 
for the e-e basis. We allow e-n indices to be zero, and 
apply a number of constraints on the linear parameters: 
(a) wc eliminate all index sets except those in which the 



k 

e-n indices are of the form M = ( q ^ 



or 



k 

1 



with 



k,l > (b) wc eliminate all parameters that do not sat- 
isfy MOD(fc - 1,3) = MOD(/ - 1,3); (c) we equate any 
two linear parameters X^^^ and A^^^ if the zeros of 
and /i^ are in the same position and their nonzero com- 
ponents satisfy INT [{ki - l)/3)] = INT [(/s2 - l)/3] and 
INT [(^1 - l)/3)] = INT [(^2 - l)/3]. These constraints 
are applied in addition to the generic symmetry and 
cusp constraints. Tabic |I] summarizes the notation for 
the terms we have introduced. 
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TABLE I. Notation for Jastrow terms correlating n electrons and m nuclei using different basis functions. 

Name Basis set Cutoff function Special constraints 

Nn^rn Natural powers Polynomial None 

Fn.m Powers of r/(r'' + a) None None 

Bn,m Powers of r-/(r + o) None Boys-Handy-style indexing 

^njin Natural powers Anisotropic polynomial None 

C„,m Cosines None None 

Ki,m Natural powers times unit vectors Polynomial Dot product 



IV. RESULTS 

In the present work we have used a variety of methods 
to optimize our Jastrow factors, namely variance mini- 
mization, minimization of the mean absolute deviation 
of the local energy with respect to the median energy, 
and linear least-squares energy minimization psl . l26l | . All 
of our final wave functions are energy-minimized except 
where otherwise stated. Starting with the Hartree-Fock 
(HF) wave function, we progressively introduce Jastrow 
terms and re-optimize all of the parameters simultane- 
ously. Optimizing the Jastrow factor term-by-term is 
unnecessary in practical applications, but here it allows 
us to understand the importance of the different terms. 
We refer to the total number of optimizable parameters 
in the wave function as Np. 

The correlation energy is defined as the difference be- 
tween the HF energy and the exact energy, £'hf — £-oxact- 
The fraction of the correlation energy retrieved in a VMC 
calculation with a given trial wave function 5*, 



/ce[*] = 



VMC 



EuF — £^cxact 



(35) 



is a measure of the quality of ^. We refer to the dif- 
ference between the DMC and HF energies as the DMC 
correlation energy, E'hf — -E-dmc [^] • The fraction of the 
DMC correlation energy retrieved in VMC, 



/dce[*] 



E 



HF 



VMC 



DMC 



(36) 



measures the quality of the Jastrow factor, since a perfect 
Jastrow factor would make the VMC and DMC energies 
coincide We define the fraction of the remaining 

DMC correlation energy recovered by a wave function 
\E'2 with respect to another as 



£^VMc[^l] — £^VMc[^2] 
i^VMc[*l] - £^DMC[*2] 



(37) 



The local energy of an electronic configuration R is de- 
fined as EiiR) = 4'T\R)i?(R)*T(R), where H{K) is 
the Hamiltonian operator. The variance of the local en- 
ergies encountered in a VMC calculation, which we shall 
refer to as the VMC variance, tends to its lower bound 
of zero as tends to an eigenstate of the Hamiltonian, 
and is thus a measure of the overall quality of the trial 
wave function. 



A. Homogeneous electron gases 

1. One- dimensional homogeneous electron gas 

We have studied a ID HEG of density parameter 
rs — 5 a.u. consisting of 19 electrons subject to peri- 
odic boundary conditions using a single Slater determi- 
nant of plane- wave orbitals. The ground-state energy of 
an infinitely thin ID HEG in which electrons interact by 
the full Coulomb potential is independent of the mag- 
netic state, and hence we have chosen all the electrons to 
have the same spin. This system is unusual in that the 
nodal surface of the trial function is exact, and therefore 
DMC gives the exact ground-state energy, which we have 
estimated to be —0.2040834(3) a.u. per electron. Excel- 
lent results were reported for this system in Refs. [28l[29| 
using wave functions with e-e backflow transformations 
[1, which preserve the (exact) nodal surface of the 
Slater determinant. 

We have investigated the improvement in VMC results 
when various terms are added to an e-e Jastrow factor 
J = A^2,o+C2,Oi both with and without backflow transfor- 
mations. In the absence of backflow. we find that includ- 
ing N^fi, Csfi, or A^3_o + C3,o improves the VMC energy, 
while the subsequent addition of C4.0 yields no further 
gain. VMC gives an almost exact energy with backflow 
and J = N2,o, and therefore no further reduction is pos- 
sible by including more Jastrow terms. However, the 
addition of N^^ + C^^ reduces the VMC variance by a 
factor of five, giving a variance that is an order of magni- 
tude smaller than that reported in Ref. [1^ for a similar 
calculation. 



2. Two-dimensional homogeneous electron gas 

We have studied a paramagnetic 2D HEG with 42 elec- 
trons per simulation cell at = 35 a.u., which lies close 
to the Wigner crystallization densitypredicted by Drum- 
mond and Needs [13 . Kwon et al. [1^ found that three- 
electron correlations are important at low densities, and 
that the effect of a three-electron Jastrow factor on the 
VMC energy is comparable to that of backflow. This 
makes low densities appealing for testing higher-rank Jas- 
trow terms. 

The VMC energy and variance obtained using different 
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I 1 I 1 I 

2x10"^ 4x10'^ 6x10 
VMC variance (a.u.) 



FIG. 1. (Color online) VMC energies Eymc against the VMC 
variance for the 2D HEG at Vs — 35 a.u. using different Jas- 
trow factors, along with the DMC energies for reference. The 
error bars are smaller than the size of the symbols, and "(BF)" 
indicates the use of backflow. 



Jastrow factors with and without backflow is plotted in 
Fig.[T] The addition of an N^^ term to J = 7V2,o recovers 
81% of the remaining DMC correlation energy without 
backflow and 49% with backflow. The €2,0 term further 
reduces both the VMC energy and variance. The use of a 
Ca^o term recovers 10% of the remaining DMC correlation 
energy when added to J = iV2,o + C2,o, but it was not 
used further since the lack of a cutoff function makes 
calculations with Ca n too costly for the little benefit it 
provides. 

We have also performed DMC calculations using two 
different Jastrow factors in the presence of backflow in 
order to quantify the indirect effect of the quality of the 
Jastrow factor on the quality of the nodes of the wave 
function. We obtain a DMC energy of -0.0277072(1) a.u. 
per electron for J = ^-nd a lower energy of 

-0.0277087(1) a.u. per electron for J iV2,o+A^3,o+C'2,o. 
This supports the idea that a better Jastrow factor al- 
lows the backflow transformation to shift its focus from 
the "bulk" of the wave function to its nodes, thus im- 
proving the DMC energy. 



radial grid. We have investigated the use of Jastrow fac- 
tors with up to four-body terms, but we have not used 
backflow for these systems. The energies of Chakravorty 
et al. HI] have been used as "exact" reference values. 

We obtain lower single-determinant VMC energies for 
the Be, B, and O atoms with J = ^2,0 + ^i,i + ^2,1 than 
reported in Refs. [U, [s^l. We obtain further small im- 
provements in the VMC energies by including cither F^^ 
or i^3_i Jastrow terms, but their combination, -Fs.o + ^a,!, 
is not found to be advantageous over using the terms in- 
dividually. This indicates that F^^q and F3.1, the latter of 
which provides a slightly lower VMC energy than the for- 
mer, have nearly the same effect in atoms. These three- 
electron terms should be particularly useful in describing 
correlations involving electrons in the atomic core region. 
We expect F^^i to be more useful than ^3^0 in molecules 
and solids because it should be able to adapt to the dif- 
ferent length scales in these systems, whereas ^3^0 offers 
a homogeneous description of three-electron correlations. 
We have investigated the effect of adding a F4.1 term in 
Be and O, but it does not reduce the VMC energy or 
variance when added to J = ^2,0 + ^"1,1 + ^2,1 + ^3,1- 

Our best VMC energies of -14.6522(1) a.u. (Be), 
-24.6309(2) a.u. (B), and -75.0381(3) a.u. (O), cor- 
respond to fractions of the DMC correlation energy of 
94.0(1)%, 91.8(1)%, and 94.6(1)%, respectively. The 
VMC energies are compared with the best single- 
determinant nonbackflow VMC values we could find in 
the literature in Table Hill 



C. BeH, N2, H2O, and H2 molecules 

The BeH, N2, H2O, and H2 molecules are strongly in- 
homogeneous and anisotropic systems. We have used ba- 
sis sets of moderate quality for the single-electron orbitals 
of BeH and N2 in order to investigate the extent to which 
the Jastrow factor can compensate for the deficiencies of 
the basis sets, especially via one-electron terms Ni^m- For 
H2O and H2 we have used very good basis sets. We have 
also tested anisotropic Jastrow factors in N2, and a van 
der Waals-like Jastrow factor for H2. 



B. Be, B, and O atoms 

While excellent descriptions of these atoms can be ob- 
tained within VMC and DMC using multideterminant 
wave functions with backflow correlations [U, [3l|, we 
have used single-determinant wave functions since we are 
only interested in the effects of the Jastrow factor. We 
have studied the ground states of the Be, B, and O atoms, 
corresponding to ^S, ^P, and electronic configurations, 
respectively. The ATSp2k code [s^l was used to gener- 
ate numerical single-electron HP orbitals tabulated on a 



1. BeH molecule 

We have studied the all-electron BeH molecule in the 
^S"'' ground state configuration at a bond length of 
2.535 a.u. [1^. We have used a single-determinant wave 
function containing Slater-type orbitals generated with 
the ADF package [30], with which we obtain a reference 
DMC energy of -15.24603(4) a.u. 

The addition of Ni^2 to J = A^2,o + -^i,i + -^2,1 recovers 
11% of the remaining DMC correlation energy. We find 
no significant gain from adding either an A^2,2 term or an 
Ns^i term to J = 7V2,o + ^1,1 + N2,i + iVi,2. 
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2. N2 molecule 



We have studied the 



ground state of the N2 
molecule at the experimental bond length of 2.074 a.u. 
(35I I HF orbitals were generated in a Slater-type basis us- 
ing the ADF package [3g ■ Our VMC results for different 
Jastrow factors are given in Table along with relevant 
reference energies. 

Adding an term to J = iV2.o + A^i.i + A^2,i factor 
recovers 33% of the remaining DMC correlation energy 
and leads to a significant reduction in the VMC variance. 
The subsequent addition of N2.2 provides a reduction in 
the VMC energy of 13% of the remaining DMC correla- 
tion energy. We have tested adding iVa.o, -/V3.1, and A'4^q 
terms to J = A^2,o + -^1.1 + ^2.1 + -^2,2, but neither of 
these yield any improvements in the VMC energy. 

The anisotropy of this system is expected to be cap- 
tured by terms containing e-n functions that treat the 
bond as a special direction. We have aligned the z-axis 
of our reference frame along the N-N bond in our calcula- 
tions, and A\ ^ is then the simplest explicitly anisotropic 
term that reflects the geometry of the system. The ^ 
and A\ ^ terms must be zero by symmetry and we have 
therefore not used them. There are five spherical har- 
monics with Z = 2, which are respectively proportional 
to xy, xz, yz, x^ — y^, and —x^ ^ 2/^ + 2z^. In our calcu- 
lations we find that only the last one of these, which we 
refer to as z^, has a significant effect on the VMC energy. 

The VMC energy with J = iV2,o + A^i,i + A^2,i + 
A\ I is within statistical uncertainty of that with J = 
-^2,0 + A^i.i + -^2,1 + -^1^2, but the former Jastrow factor 
contains about a third fewer parameters than the lat- 
ter. The combination of the A^i,2 and A\ ^ terms into 
J = A''2,o + A^i.i + N2.1 + A^i,2 +^11 does not improve 
the VMC energy compared with the other two Jastrow 
factors. These results suggest that the terms A^i^2 and 
A\ I play similar roles in the wave function, which we 
find reasonable since iVi.2, although constructed from 
isotropic basis functions, contains the right variables to 
capture the symmetry of the molecule in much the same 
way as ^ does. We have plotted the ^ term for 
J = N2fi + Nis + N2,i + Al^^ and the Ni^2 term for 
J = N2fi + Ni^i + N2,i + iVi,2 'm Fig. H where the sim- 
ilarity between the terms can be seen. The value of the 
iVi.2 term is roughly the same as that of Af ^ offset by 
a positive amount, and this shift is likely to be compen- 
sated for by the other Jastrow factor terms. Both terms 
increase the value of the wave function in the outer region 
of the molecule with respect to that in the bond region. 

We have added different combinations of anisotropic 
terms to J = A^2.o + -^1.1 + -^2.1- The e-e-n A2 1 term 
retrieves less correlation energy than the e-n j term. 
Combining terms with spherical harmonics of / = 1 and 
I = 2 improves the VMC energy significantly with respect 
to using I = 1 only. The anisotropic Jastrow factor J = 
N2,o + Ni^i + N2.1 + ^f,i + Aj'i + + A(\, which 
contains up to e-e-n correlations and has 191 optimizable 
parameters, recovers 93.3(1)% of the DMC correlation 



energy, which is greater than the 93.0(1)% retrieved by 
our best isotropic Jastrow factor J = iV2,o + -^i,i + -^2,1 + 
Ni,2 + ^^2,2 + A^s.o, which includes more costly e-e-n- 
n and e-e-e correlations and contains 260 optimizable 
parameters. We conclude that anisotropic functions are 
an important tool in the construction of compact Jastrow 
factors for strongly anisotropic systems. 

Toulouse and Umrigar obtained 90% of the DMC cor- 
relation energy with a single-determinant wave function 



in Ref. [3J] , and with our best Jastrow factor we retrieve 



93% of the DMC correlation energy. We have also opti- 
mized a single-determinant backflow wave function with 
our best Jastrow factor and we obtain a VMC energy of 
— 109.4820(6) a.u. (89% of the correlation energy), which 
is of similar accuracy to the multidetcrminant VMC en- 
ergy of -109.4851(3) a.u. (89.6% of the correlation en- 
ergy) obtained by Toulouse and Umrigar. 



3. H2 O molecule 

Single-particle spin-unrestricted HF orbitals for the 
^Ai ground state of H2O were generated using the CRYS- 
TAL Gaussian basis set code [Sal- The basis set for O 
contains 14 s-, 9 p-, and 4 d-orbitals, and that for H con- 
tains 8 S-, 4 p-, and 3 d-orbitals. Electron-nucleus cusps 
have been added using the scheme of Ma et al. [s^ . We 
have simulated a water molecule with a bond length of 
rpH = 1.8088 a.u. and a bond angle of Zhoh = 104.52° 

Adding an 7Vi_2 term to J = N2fi + Ni^i -\- A^2,i gives 
only a very small improvement for H2O, compared with 
the more substantial improvements obtained with this 
e-n-n term for BcH and N2. The iVi.2 term acts as a 
correction to the single-electron orbitals, and we believe 
that it is unimportant in this case because wc have used 
very accurate HF orbitals, whereas the single-electron 
orbitals used for BeH and N2 are considerably less ac- 
curate. We find additional small improvements to the 
energy of H2O from adding N^^q and N^^i terms to 
J = N2fi + Ni^i+ N2,i. 

Clark et al. obtained 92% of the DMC correlation en- 
ergy with a single-determinant wave function in Ref. [4l| , 
and with our best Jastrow factor we recover 95.5% of the 
DMC correlation energy. 



4-. H2 molecule 

The energy of the first triplet excited state ('^S^) of 
H2 has a very shallow minimum corresponding to a large 
bond length of nearly 8 a.u. Although the exchange 
interaction falls exponentially with increasing internu- 
clear separation, Kolos and Wolniewicz found that it con- 
tributed significantly to the energy even at the large dis- 
tance of 10 a.u. [42] The strong interplay between the 
attractive dispersion forces and the repulsive exchange 
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TABLE II. VMC energies (E) and variances (V) for the N2 molecule using diflterent Jastrow factors, including explicitly 
anisotropic terms. We have used a bond length of tnn = 2.074 a.u [35l ]. 



Np E (a.u.) V (a.u.) /ce (%) /dce (%) 



jir iimiL iroiQ rv,er. 




1 Oft QQ9Q 

— iUo.yyzy 




n 
U 


u 




18 


— 109.102(1) 


5.275(4) 


19.9(2) 


21.3(2 


A^2,o+A''i,i 


27 


— 109.3739(6) 


3.681(3) 


69.4(1) 


74.3(2 


TV 7" 1 A r 1 TV 7 

iV2,0+A^l,l+^1.2 


49 


— 109.3796(6) 


3.595(2) 


70.4(1) 


75.4(2 


TV T' 1 A r 1 71 T" 

N2,0+NiA+N2,l 


80 


— 109.4441(4) 


1.667(2) 


82.16(7) 


87.9(1 


Af^ r. -L- Af-, -, -1- Af^ -, -1- A/"-, -1 


1 09 
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A''2 0+A^1 1 +A^2 1 +A^1 2+A^2 2 
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-109.4697(4) 


1.088(3) 


86.82(7) 


92.9(1 


A''2 , + A^i , 1 + iV2 , 1 + iVi , 2 + A''2 , 2 + iVa , 


260 


-109.4702(3) 


1.083(2) 


86.91(5) 


93.0(1 




36 


—109.3770(6) 


3.670(2) 


69.9(1) 


74.9(2 




89 


-109.4660(3) 


1.116(2) 


86.14(5) 


92.2(1 


A^2,0+A^l,l+iV2,l+^l,l+Af'i 


97 


-109.4669(3) 


1.073(2) 


86.31(5) 


92.4(1 


iV2,o+Afi,i+A^2,i+^!,i+Ai,i 
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-109.4707(3) 


1.072(2) 


87.00(5) 


93.1(1 


iV2,()+Afl,l+A^2,l+^l,l+^l'l+-42,l+^ii'l 


191 


-109.4714(3) 


1.036(4) 


87.13(5) 


93.3(1 


SD VMC from Ref. [34|'' 




-109.4520(5) 




83.59(9) 


89.5(2 
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-109.5060(7) 




93.4(1) 


100.0(2 


Exact from Ref. [37] 




-109.5421 
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^ For rjMN = 2.075 a.u. We do not expect that this small difference in bond length will aflect the 
comparison between energies significantly. 




FIG. 2. (Color online) Plots of the Af^i term (left) and A''i,2 term (right) of different Jastrow factors for N2 as a function of 
the position of an electron in a plane containing the nuclei (black circles). 



interaction requires that both be accounted for to af- 
ford an accurate description of the triplet state. This 
makes the system appealing for studying the construc- 
tion of four-body Jastrow factor terms to describe van 
der Waals-like interactions. 

We generated numerical HF orbitals tabulated on an 
elliptical grid using the 2dhf package [i^. HF theory 
predicts no binding for the triplet state at any separa- 
tion, and therefore any binding that occurs in VMC can 
be attributed to the Jastrow factor. We have studied 
the H2 molecule in its triplet spin state at the internu- 
clear distance of 7.8358 a.u. This separation and the 
corresponding energy of —1.0000208957 a.u. were found 
by fitting a quadratic function to the data of Staszewska 
and Wolniewicz [i^ . 



Previous QMC calculations on H2 at different inter- 
atomic distances have used Jastrow factors with up to 
four-body correlations where the cusp conditions were 
not enforced [3, [13 : instead relying on the variance min- 
imization method to find parameter values that approx- 
imately satisfy the cusp conditions. This was found to 
be advantageous for this system because the additional 
variational freedom yielded a better description in VMC 
than when the cusp conditions were obeyed exactly [l^] ■ 
The violation of the cusp conditions is potentially catas- 
trophic in DMC calculations, but previous studies have 
restricted the use of such terms to VMC. 

For H2 we have optimized Jastrow factors consisting 
of the single e-e-n-n terms 1/2,2, ^ll^i ^^'^ ^2.2 (see Ta- 
ble m at several expansion orders, where no constraints 
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are enforced at e-e or e-n coalescence points. We have 
used variance minimization for these Jastrow factors as 
we found that it produces better results than energy min- 
imization. We have also optimized Jastrow factors con- 
sisting of different sums of terms which satisfy the cusp 
conditions using energy minimization. The results are 
shown in Fig. |31 
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FIG. 3. (Color online) Difference between the VMC and ex- 
act energy against the number of wave function parameters 
for the H2 triplet ground-state using different Jastrow fac- 
tors. Only the "multiterm" Jastrow factor enforces the cusp 
conditions. The error bars are smaller than the size of the 
symbol where not shown. All of the energies in this plot are 
lower than —1 a.u. and therefore the wave functions predict 
binding. 



the DMC correlation energy, offering a good description 
of the system without reaching the accuracy of the more 
generic ^ and i?2,2 terms. 

A V2,2 term without e-e functions consists of contribu- 
tions proportional to rn ■ Vjj, where the prcfactors de- 
pend explicitly on rn and r^j, and implicitly on r/,/, and 
this functional form is that of a dipole-dipolc interaction 
potential. Our best such V2.2 term retrieves 69% of the 
DMC correlation energy, which amounts to 0.0000262(3) 
a.u., and we regard this as a measure of the pure van der 
Waals correlation energy of this system. 

The multiterm Jastrow factors contain the usual A^2,0j 
-^1,2, and N2.1 terms, and for each combination of 
these we have added a V2,2 term without e-e functions 
obeying the cusp conditions to study its effect. J = iV2.o 
retrieves 44% of the DMC correlation energy, and adding 
the V2,2 term retrieves 85% of the remaining DMC cor- 
relation energy. The effectiveness of ¥2^2 progressively 
drops as more terms are added, and it retrieves 43% of 
the remaining DMC correlation energy when added to 
J = N2fi + Ni^i + A^2,i + -^1,2- In all cases, V2,2 is found 
to lower the VMC energy by a larger amount than any 
of the Nn^m terms. 

Our best multiterm cusp-enforcing Jastrow factor re- 
trieves 97% of the DMC correlation energy with 77 wave- 
ftmction parameters, comparable with the 98% retrieved 
with the cusp-violating i^^l^ ^^'^ ^2,2 terms with a sim- 
ilar number of parameters. For larger systems where van 
der Waals interactions are important, we expect the vi- 
olation of cusp conditions to cause statistical problems, 
and the ^2,2 term would become an effective way of im- 
proving the description of the system in a multiterm Jas- 
trow factor. 



We have performed the DMC calculations using our 
best i?2,2 Jastrow factor and obtain a reference DMC 
energy of —1.0000207(1) a.u. We have not encountered 
any statistical problems in the DMC calculations with 
this cusp-violating wave function. Such issues can occur 
when the local energy presents a negative divergence in 
a region of configuration space with a significant proba- 
bility of being sampled. We have verified that our wave 
function causes a negative divergence in the local energy 
when an electron coalesces with a nucleus, a point at 
which the wave function is likely to be relatively large. 
We therefore conclude that the region of influence of this 
divergence is sufficiently small that statistical problems 
do not arise in practice. 

The F2^^ and -82,2 terms only differ in that the lat- 
ter uses Boys-Handy-style indexing, which yields slightly 
lower VMC energies than standard indexing in most cases 
for a fixed number of parameters. Our best -Ff'f ^ and 
B2,2 Jastrow factors retrieve 99% of the DMC correla- 
tion energy in VMC. 

The V2,2 term is designed to describe van der Waals 
correlations, and contains e-e functions which introduce 
other correlations. Our best V2,2 term recovers 92% of 



D. Discussion of molecular results 

In Fig. m we have plotted the fraction of the DMC 
correlation energy retrieved by different Jastrow factor 
terms for the BeH, N2, H2O, and H2. Our purpose is 
to visualize the importance of different terms in different 
system, and to this end we do not include anisotropic or 
cusp-violating terms, and we only consider the addition 
of terms in a specific order. 

The N2fi term represents the simplest description of 
electronic correlations and typically retrieves 20-25% of 
the DMC correlation energy. This e-e term greatly dis- 
torts the charge density of the HF wave function, and 
the iVi.i term repairs this, retrieving an additional 45- 
50% of the DMC correlation energy. In the case of the 
more diffuse II2 molecule the A'2,0 and A^i^i terms have 
a different relative importance, but J = N2.0 + A^i.i re- 
trieves 70-75% of the DMC correlation energy in the four 
molecules. 

Like iVi^i, Ni^2 acts as a correction to the single- 
electron orbitals. This term provides no significant bene- 
fit in H2O, where we have used high-quality orbitals, but 
it recovers 7% of the DMC correlation energy in II2 . 
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FIG. 4. (Color online) Fraction of the DMC correlation en- 
ergy retrieved by different Jastrow factor terms, added in the 
specific order depicted in the graph (starting from the bot- 
tom), for the molecules studied in this paper. 



The effect of Ni^2 in N2 is noteworthy in that the 
energy reduction obtained by adding this term to J = 
-^2,0 + -^1,1 is about a factor of four times smaller than 
when added to the more accurate J = ^^2. + ^1.1 +^2.1- 
One would expect a term to retrieve more correlation en- 
ergy when added to a smaller Jastrow factor, and this is 
the case for Ni^2 in the other molecules. We think that 
the distortion in the charge density caused by N2,i in 
N2 is such that the single-electron correction effected by 
.^1,2 becomes more useful in its presence. 

The A^2,i term added to J = iV2,o+.^i,i+^i,2 captures 
an additional 15-20% of the DMC correlation energy. 
Higher-order terms added to J = iV2,o + -^1,1 + .^2,1 + 
■^1,2 yield significant gains in relative terms, with e-e-n- 
n terms retrieving 13% and 43% of the remaining DMC 
correlation energy remaining for N2 and H2 , respectively, 
and the e-e-e-n term recovering 17% of the remaining 
DMC correlation energy. 



tors. We have also introduced anisotropic cutoff func- 
tions. The formalism may be applied to systems with 
particle types and external potentials other than elec- 
trons and Coulomb potentials. 

Optimization of the wave function is one of the most 
human- and computer-time consuming tasks in perform- 
ing QMC calculations. We have performed tcrm-by-term 
optimizations to understand how different terms in the 
Jastrow factor contribute to the electronic description of 
a system, and we hope that our analysis will serve as a 
guideline for constructing Jastrow factors for other sys- 
tems. 

We have tested these terms on HEGs, atoms, and 
molecules. The variational freedom from the higher-order 
terms generally improves the quality of the wave func- 
tion. We have only considered single-determinant wave 
functions in this study, although our Jastrow factor can 
of course be used with other wave function forms. 

We have demonstrated the construction and applica- 
tion of an e-e-n-n Jastrow factor term designed to de- 
scribe van der Waals interactions between atoms. This 
term retrieves a large fraction of the van der Waals cor- 
relation energy in tests on the triplet state of H2 at the 
proton separation that minimizes the total energy of the 
system. 

We have found evidence for the importance of three- 
electron Jastrow terms in the low-density ID and 
2D HEGs. Improving the Jastrow factor for single- 
determinant backflow wave functions also leads to small 
improvements in the DMC energy of the 2D HEG. This 
demonstrates the indirect effect that improving the Jas- 
trow factor can have on improving the nodal surface, as 
reported in Ref. [l3j . 

We have made efforts to obtain accurate single- 
determinant VMC energies for most of the systems stud- 
ied, but for BeH and N2 we deliberately used inferior 
one-electron basis sets to see whether we could compen- 
sate for this with one-electron Jastrow terms. We find 
that this goal can be achieved by including an Ni^2 Jas- 
trow term or anisotropic c-n terms, along with the usual 
A^i,i term. 



E. Summary of results 

Table lllll gives a comparison of the best single- 
determinant nonbackflow VMC energies we have found 
in the literature with those obtained in this work. 



V. CONCLUSIONS 

We have described a generalized Jastrow factor al- 
lowing terms that explicitly correlate the motions of n 
electrons with m static nuclei. These terms can be 
parameterized using various basis sets, including terms 
that involve dot products of interparticle position vcc- 
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Appendix A: Implementation 

In this section we describe our design choices in imple- 
menting our Jastrow factor in the CASINO code [HI . The 
implementation principles are modularity and extensibil- 
ity, embracing the flexibility that this Jastrow factor has 
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TABLE III. Best single-determinant nonbackflow VMC energies (a.u.) found in the literature and those from this work, along 
with single-determinant DMC and exact energies for reference. 



System 



This work 



Literature 



DMC 



Exact 



ID HEG (r«=5, iV=19) 

2D HEG {rs=35, iV=42) 

Be 

B 

O 

BeH 

N2 

H2O 

H2 



-0.2040833(2) 
-0.0276112(6) 

-14.6522(1) 

-24.6309(2) 

-75.0381(3) 

-15.2412(3) 

-109.4714(3) 

-76.4068(2) 
-1.00002045(3) 



-14.64972(5)" 
-24.62936(5)" 
-75.0352(1)" 
-15.228(1)*^ 
-109.4520(5)^= 
-76.3938(4)'^ 



-0.2040834(3) 
-0.0277087(1) 
-14.65717(4) 
-24.64002(6) 
-75.0511(1) 
-15.24603(4) 
-109.5060(7) 
-76.4226(1) 
-1.0000207(1) 



-0.2040834(3) 

-14.66736 
-24.65391 
-75.0673 
-15.2482 
-109.5421 
-76.438 
-1.0000208957 



"Ref. 
•^Ref. 

Ref. 
•^Ref. 



34| 
45| 
341 
[411 



(using a slightly different bond length). 



by design rather than impeding it by focusing too strictly 
on performance. 



1. Basis functions 

The most important design requirement for modularity 
is that basis sets be dealt with separately rather than in- 
cluded in the Jastrow factor code. Basis functions, along 
with their derivatives when required, are evaluated and 
stored for later use. For any given term with expansion 
orders p and q there are pN{N — l)/2 e-e and qNM 
e-n functions to evaluate and store. Single-electron up- 
dates involve recalculating p{N — 1) and qM of these 
functions. Note that the number of basis functions that 
must be evaluated is independent of the term ranks n and 
TO. Furthermore, it is possible to allow different terms 
to share basis functions that do not contain optimizable 
parameters so that, e.g., the natural powers involved in 
computing A^2,o and iVi^ can be re-used for N2s- 

A number of properties of the basis functions arc re- 
quired to construct the Jastrow factor. Equation re- 
quires knowledge of whether basis functions are symmet- 
ric or antisymmetric, and their value, first radial deriva- 
tive, and angular dependency at the origin are required 
by Eqs. (Uni), (EH), (112]), dMl), dH, and i^. The one- 
contribution, no-contribution, and equal-product con- 
straints in these equations require a table indexing dis- 
tinct products of two basis functions at any value of their 
arguments. We implement interfaces that make these 
properties available so that basis functions can be treated 
as abstract objects in the construction of the Jastrow fac- 
tor, which makes implementing new basis sets straight- 
forward. 

Cutoff functions are dealt with as additional basis sets 
with an expansion order of one, and we store informa- 
tion identifying the cutoff functions that are strictly zero, 
which we use to speed up evaluation of the Jastrow fac- 



tor. 



2. Evaluation of the Jastrow factor 

For the evaluation of an arbitrary-rank Jastrow factor 
term it is necessary to use efficient procedures that iterate 
from a given set of electronic and nuclear indices i and I 
to the next in a specific order; explicit loops over scalar 
integer indices are not an option in static code since the 
loop depth is variable, and the memory usage of precom- 
puting all possible i and I scales badly with system size 
for high ranks n and to. These procedures should take 
into account which cutoff functions are zero so that par- 
ticle sets that do not contribute to the Jastrow factor are 
skipped. Efficient handling of localized Jastrow factor 
terms is important because it allows the cost of evalu- 
ating a term of any rank to scale linearly with system 
size if the cutoff lengths are held fixed. We implement 
a scheme where we construct a list of the electrons that 
are "connected" to each electron and each nucleus via 
nonzero cutoff functions. For the electronic indices, the 
value of ii is iterated between 1 and N , then the value 
of 12 is iterated over the values in the list associated with 
electron ii that are greater than ii, then the value of 13 
is iterated over the values in the intersection between the 
lists associated with ii and 12 that are greater than 12, 
and so on. The procedure that iterates over nuclear in- 
dices selects sets of nuclei whose connected-electron lists 
have nonzero intersections. We iterate over I in the out- 
ermost loop so that we can feed the intersection of all 
e-n lists into the i iterator as an initial list for index ii. 

The signature {P(i), S(i, I)} of each group of particles 
is computed inside the electronic and nuclear loops, iden- 
tifying the linear parameter channel associated with the 
group of particles. Wc then loop over linear parameters 
in the channel, computing the products of the relevant 
basis functions which have already been precomputcd. 
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In terms without indexing constraints, consecutive linear 
parameters tend to have very similar expansion indices, 
i.e., they multiply most of the same basis functions. In 
order to save multiplications, it is convenient to buffer 
partial products so that, e.g., if only the last two of six 
expansion indices change from one parameter to the next, 
we can recover the product of the first four functions from 
the previous index set and save three of the five multi- 
plications required to combine the six basis functions. 

In typical QMC calculations individual electron moves, 
rather than full configuration moves, are proposed, which 
requires computing an acceptance probability involving 
the ratio of the trial wave function at the proposed and 
original positions. To calculate this efficiently one needs 
to compute the part of the Jastrow factor which depends 
on the position of a single particle i, ignoring the contri- 
butions not involving i. In our implementation we evalu- 
ate this one-electron Jastrow factor using Eq. ([S]), where 
we fix ii = i and iterate over the rest of i. The main 
difference from the evaluation of the full Jastrow factor 
is that i is not sorted, and the permutation required to 
sort i, which amounts to inserting i at the correct posi- 
tion in i, needs to be taken into account in the presence 
of antisymmetric e-e basis functions, since sign changes 
will be required by Eq. ([T^ . The evaluation of gradi- 
ents and Laplacians of the Jastrow factor term can be 
easily accommodated in the one-electron Jastrow evalu- 
ation code. 

For performance reasons it is advisable to implement 
versions of the evaluation procedure for fixed ranks, with 
fixed-depth loops that can be optimized by compilers. 
We implement three optimized versions: one for e-e 
terms, one for e-n terms, and one for e-e-n terms. Other 
terms are handled by three generic procedures: one for 
terms without e-e functions, one for terms without e- 
n functions, and one for terms with both e-e and e-n 
functions. 



Appendix B: Term construction example 

Let us consider the iVa o term used for the 2D HEG in 
Section flV A 21 This is a system with iV = 42 electrons, 
half of each spin, and the A^a.o term corresponds to n = 3, 
771 = 0, expansion order p = 4, and basis functions 



Lf 



Li 



(Bl) 



The spin-pair dependency matrix P is of size N x N, 
but in practice we specify a reduced 2x2 version of this 
matrix where each row (column) corresponds to different 
spins. For this system this matrix is 



P = 




(B2) 



where we assign parallel-spin electron pairs a spin-pair 
dependency index of 1, and antiparallel-spin electron 



pairs an index of 2. The P matrix can be regarded as 
an intrinsic property of the system, but in some cases 
additional symmetries can be imposed in order to reduce 
the number of parameters in the Jastrow factor term; for 
example we would achieve this by setting all elements of 
P to 1 in the 2D HEG, choosing to ignore the distinction 
between parallel- and antiparallel-spin electron pairs. 

The elements in P determine pairwise properties and 
objects. For example, there are as many sets of nonlinear 
parameters in a Jastrow factor term as different values 
in P; in this case, there are two cutoff lengths, Li for 
parallel-spin electron pairs and L2 for antiparallel-spin 
electron pairs. 

Four distinct types of three-electron groups can be 
formed: three up-spin electrons (ttt)j two up-spins and 
one down-spin (tT-l)! up-spin and two down-spins 
(til), and three down-spin electrons (4-4-4.). The spin- 
pair dependency matrices for these four groups are 




p(nt) 



1 2 
p(m) = 1 1 2 
,2 2 0, 




(B3) 



'0 2 2' 
; P(m) = I 2 1 I . (B4) 
,2 1 0, 



The matrices P(ttt) a-nd P(i4'4') a-re identical, and there- 
fore correspond to a linear parameter channel which is 
used for groups of three electrons with parallel spins. The 
matrix P(t4'4') can be transformed into ^.(tti) via, e.g., 
the permutation matrix 



U: 




(B5) 



and therefore both matrices correspond to a second lin- 
ear parameter channel which is used for groups of three 
electrons with mixed spins. The signature of the first 
channel, which we refer to as the ftt channel, is P(ttt): 
and the signature of the tt4' channel is P(tti)- Both 
of these matrices are considered sorted by our matrix- 
sorting algorithm. 

The symmetry constraints for the parameters in each of 
the two channels depend on the above matrices. P(ttt) 
is invariant with respect to the application of any per- 
mutation, and for the ttt channel Eq. equates any 
two parameters with the same indices in a different or- 
der. The signature of the tt4' channel is only invariant 
with respect to one nontrivial permutation. 



U: 




(B6) 
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and therefore the symmetry constraints for this channel 
are 



A 



/O 1 2\ 
1 2 
v2 2 Oy 
' Vij 




) 



= A 



/O 1 2\ 
1 2 
V2 2 0/ 
^ i^ij 




; 



(B7) 



We do not aUow expansion indices to be zero in the 
iVs^o term, and the presence of cutoff functions prevents 
the "one contribution" constraint of Eq. (|20p from being 
satisfied. It is therefore not possible to impose the Kato 
cusp conditions on the A^3,o term. This is not a problem 
in practice because we use this term in conjunction with 

Since the basis functions are isotropic, the only con- 
straint applicable to this term at coalescence points is 
Eq. (|22p for Z = m = 0. The derivative of the radial 



projection of the basis function is 



dr 



C 



(B8) 



J r=0 



and the constraint equation for index Vij is 



o.p. 

E 



/ P^k\ 



R, P, 



\P^k Pjk / 



C 



/ Vij i^tk] ^P.j 



1/ 



'jk 



\Vik i^jk / 



= 



(B9) 

Any two products of pairs of natural powers is equal if the 
sum of the exponents in each of them is equal. Therefore 
the "equal-product" constraint for the N^^q term is 



o.p. min(/— l.p) 

E = E 



(BIO) 



where I ranges from 2 to 2p, and the constraint equation 
for index Vi^ thus reduces to 



min(Z — l.p) 

E 



( ° 

\p,k 

( 
2 



P^J P^k\ 
P,k 

P,k J 

2 i^ik ^ 
I - Uik 
\uik. l-Vtk / 



-A 



( ° 
P 

\P^k 
( 
1 



P^J P^k\ 

P,k 

Pjk / 

1 l^ik \ ^p, 

I - Vik 



C 



\iyik I - Vtk / 



(Bll) 



The constraints for the three expansion indices are equal by symmetry in the f'f'f channel, and there are two sets 

of constraints in the tt4- channel. 
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